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Color Glass Condensate (CGC) provides a classical description of dense gluon matter at high 
energies. Using the McLerran-Venugopalan (MV) model we calculate the initial energy density e(r) 
in the early stage of the relativistic nucleus-nucleus collision. Our analytical formula reproduces the 
quantitative results from lattice discretized simulations and leads to an estimate e(r = 0.1 fm) = 
40 ~ 50 GeV-fm" 3 in the Au-Au collision at RHIC energy. We then formulate instability with 
respect to soft fluctuations that violate boost invariance inherent in hard CGC backgrounds. We 
find unstable modes arising, which is attributed to ensemble average over the initial CGC fields. 



In the relativistic heavy-ion collision Color Glass Con- 
densate (CGC) describes the initial state of energetic 
gluon matter with the transverse momentum p t up to 
the saturation scale Q s which universally character- 
izes the hadron or nucleus wavefunction in the small- 
er regime 0, 0, |, i, i, H, 0, I, I, O- Given the 
scale Q s at a certain value of Bjorken's x, the gluon 
distribution probed by processes with Q 2 <C Q 2 is so 
dense that coherent fields should be more relevant than 
the individual particle picture during r < Qj 1 . Physi- 
cally Q 2 S corresponds to the transverse density of partons 
and is estimated by the Golec-Biernat and Wiisthoff fit, 
Q 2 = Q^Xq/x^A 1 / 3 , where A is the atomic number. 
We can expect Q s around 1 ~ 2 GeV for the Relativistic 
Heavy Ion Collider (RHIC) and 2 - 3 GeV for the Large 
Hadron Collider (LHC) in case of A = 197 (Au-Au col- 
lision) assuming relevant pt is ~ 1 GeV. This transient 
but still coherentgluon matter, which is often referred 
to as "Glasma" 7], should melt toward a quark-gluon 
plasma. 

The physical property of Glasma has been mainly an- 
alyzed by numerical simulations in the lattice discretized 
formulation [J, H, H, 0, B Hi • I n this paper we aim to ap- 
proach Glasma in an analytical way along a similar line 
to the near-field expansion proposed by Fries-Kapusta- 
Li [lo| . The analytical method is desirable for a deeper 
insight into the Glasma, which presumably exists up to 
T ^ Qj 1 ~ 0-1 f m m the Au-Au (central) collision at 
RHIC energy, y/s = 200 GeV/nucleon, or even longer 
depending on the interpretation of the Glasma. In par- 
ticular, the problem of early thermalization still has in- 
teresting unanswered questions [111 ]. We will specifically 
address the following; is there any unstable mode grow- 
ing around the initial CGC fields right after the collision? 
If any, it could speed up thermalization (or isotropiza- 
tion) even in the classical regime (t < Or 1 ), besides 
non-Abelian plasma instabilities [12, |Tj, [l4|, [la] which 
take place at later times. The pioneering numerical sim- 
ulation 0] suggests the existence of "Glasma instability" , 
though the literal time scale of instability seems to be 
greater than Qj 1 by three order of magnitude, probably 
because of the choice of tiny instability seeds. The delay 
in the instability onset has also been pointed out in the 
Hard Expanding Loop (HEL) approach to non-Abelian 
plasma instabilities [151 ] . 



We shall start with the boost-invariant CGC solu- 
tion and estimate an initial energy density. It is con- 
venient to adopt the Bjorken coordinates spanned by the 
proper time r = \/t 2 — z 2 and the space-time rapidity 
f] = | ln[(< + z)/(t — z)\. The radial gauge A T = is un- 
derstood throughout this work. The canonical momenta 
(chromo-electric fields) are read in this gauge as 

E i = Td T Ai , E 71 = T^drAr, . (1) 

It should be mentioned that the metric is g TT = 1, g m — 



and 
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-1 in accord with the convention 



in Refs. [9j, [l6| . The equations of motion derived from 
Hamilton's equations lead us to 

d r E l = r" 1 A^i + rDjFji , d T Ei = t~ 1 DjFjr! . (2) 

These are the basic equations for the classical description 
valid in the early stage right after the collision. The ini- 
tial condition is uniquely determined by boundary match- 
ing at singularities of the color sources p^ (x±)8(x~ ) and 
p( 2 \x±)5(x + ) representing the propagation of Lorentz 



contracted nuclei 0, , as follows 
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where and a\ are the gauge fields at r < associ- 
ated with the right-moving nucleus along the x + axis and 
the left-moving nucleus along the x~ axis [17] . It takes 
a pure-gauge form, c*i(xx) = -{l/ig)V{x 1 _)d l V^ '(x±), 
with the Wilson line defined by 



V\x ± ) = Pexp 
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dz 



9 
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for . The Wilson line for otf^ is given by replacement 
of x~ and p^'(xx) by x + and p^(xj_) in the above 
expression. We can compute the expectation value of 
physical observables by means of the average over the 
random color distribution inside nuclei using 

(pi m Hx ± )pi n) (y ± ))=g 2 p 2 S mn 5 ab 5^(x ± -y ± ) . (5) 

Here p is the only dimensionful scale in the McLerran- 
Venugopalan (MV) model and related to the saturation 
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scale Q s . We will later present all dimensionful quantities where 
in unit of /i. 

Let us evaluate the initial energy density of the 
fields ([3} at r = with the color source average ((5]). 
To do this, we need to take an average of four Wilson 
lines ~ {Vlx^V^y^Vlu^VHvxy). We can find an 
algebraic technique in the appendix of Ref. [l8| and it is 
even possible to write a formal expression down for more 
generic color structure [n|. Alter all, it turns out that 
the transverse fields are vanishing and the longitudinal 
chromo-magnetic fields, BV Q , — F 12 (q), are § 
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The number of color is N c = 3 in QCD. We defined a 
resulting from the two-point function in terms of Ci- . 
In order to make a direct comparison to the numerical 
simulation transparently, we shall make use of the lattice 
regularization, which gives 

L/2a 



CT= 972 E 



2L 2 ^— ' 2 — cos(27rn 1 a/L) — cos(2tt n 2 a / L) 

rii — l — L/2a 
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Here L is the size of the system fixed by L 2 — ttR\, 
and a is the lattice spacing. We got rid of the zero- 
mode ni = ri2 = because of global neutrality. We nu- 
merically checked that the above logarithmic form with 
adjusted by a constant c ~ 1.36 is a quite good ap- 
proximation. Some further calculations end up with 
the same amount of the chromo-electric field squared; 
(2tr(£ ( '' 0) ) 2 ) = (2tr(B^ 0) ) 2 ). As a result, we can estimate 
the initial energy density as 
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with N c = 3 substituted. This a and L dependent re- 
sult should be interpreted carefully, while the quantita- 
tive output somehow agrees with the latest simulation by 
Lappi; our estimate by Eqs. and (jHJ) yields 0.81 and 
0.90 for L/a = 500 and 700 which are close to 0.76 and 
0.88 reported in Ref. [8j. The logarithmic singularity has 
been found also in Refs. [1, [l(|. The singularity arises 
from the approximations that we regarded the colliding 
nuclei as infinitely thin in the longitudinal direction and 
that the random color distribution is uncorrelated at ar- 
bitrary microscopic scale in transverse space. 

Next, we will step away from singularity located at 
t = by the near-field expansion in terms of r, i.e., O = 



O 



(0) 
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The first-order corrections 



are vanishing, and the second-order fields are 
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which physically represent the initial longitudinal and 
second-order transverse chromo-magnetic fields. 

Using these expressions we calculate the contributions 
to the energy density of order r 2 to find the same 
amount of chromo-magnetic and chromo-electric fields 
again. Since the initial state has non-zero longitudinal 
fields, it follows that the cross terms between the zeroth 
and second-order terms give 
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where we defined 
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We dropped terms proportional to a 3 not containing x 
because x>ff when a is small. In the same approxima- 
tion the transverse fields of order t 4 (that is, r 2 -order in 
the energy density) result in 



2tr ( £, (2) £ '(2), 



( 9 2 ^ \ 2i < B h) B h). 



= -N c (N 2 -l)a-x + 0(a 3 ) . 

After all, we get the expanded series, 
-,2 „2 
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It is obvious from Eq. ([14]) that the r expansion behaves 
badly for small value of a, which is also clear by the 
dotted curve in Fig. [1] that plots Eq. (fT4"|) . 

The naive r expansion is, in fact, ill-defined. It is be- 
cause, as pointed out in Ref. 0], the energy density be- 
haves as ~ (lnr) 2 near t = when the colliding nuclei 
are infinitely thin. Therefore, the naive Taylor expansion 
around r = is meaningless. Nevertheless, we stress that 
we can derive meaningful information from Eq. fj 14[) ; we 
know that the asymptotic form ~ (lnr) 2 in the a — > 
limit and we also know that the regularized expansion 
~ c\ ln(L/<x)[ln(L/a) + C2(r/a) 2 ] with a kept finite. The 
simplest analytical function satisfying these two require- 
ments is ~ ci{ln[L 2 /(a 2 — C2T 2 )]} 2 , that means, 



i(2) 



- 2 E h) - 4 D j(.0)Fji(0) 



= ±EV 

~ 2-^(0) 



^(2) - k D m F - 



(9) 



9 



(5 2 M) 4 




\g 2 ^Lf 



(g 2 ^a) 2 + n(g 2 ^T) 2 



(15) 



3 



■3S 

N 



0.5 



U — a=0 limit 






-a 
■ 




- 






- 

Log— Ansatz 


■ 2, 2 

■ -via 


-1/7 





FIG. 1: Comparison of the energy density in case of L/a = 
700; the data with error bar is taken from Ref. [8| . The dotted 
and solid curves represent the naive expansion in Eq. (|14[) and 
the log-ansatz in Eq. (|15[) , respectively. The dashed curve 
scales as 1/r whose starting point is chosen at g 2 fir — 0.5, 
meaning that the "formation time" [|| being g 2 ^To ~ 0.5. 



inferred from Fig. [TJ This factor might depend on Aq CD , 
and thus, the numbers listed in the second and third rows 
should be considered as the upper and lower bounds. 

It is a natural choice to take the confining scale as the 
nucleon size ~ 1 fm, and the estimate of the initial energy 
density is then e(r = 0.1 fm) = 40 - 50 GeV- fm" 3 . This 
value is significantly smaller than the previous estimates, 
130 GeV- fm" 3 in Ref. [| and 260 GeV- fm" 3 in Ref. [nj, 
reflecting difference between the choices Aq CD = 1 fm 

and Aqq D = L — 12 fm, but rather consistent with the 
simulation with color neutrality of finite nuclei taken into 
account 5] that found e(r = r D ~ 0.3 fm) = 7.1-40 GeV 
fm" 3 . 

When g 1 [it becomes larger, the energy density comes 
to scale as ~ t /t because of (almost free streaming) 
longitudinal expansion 0, [ll[ . It should be noted that 
the scaling law in the classical regime is different from the 
(one-dimensional) hydrodynamic one — (to/t) 4 / 3 . For 
reference we plot the scaling behavior e(T)/e(ro) = tq/t 
in Fig. [T] indicated by the dashed curve with a choice 
of g 2 [J,To = 0.5. The expanding system at late times is 
dilute so that this scaling behavior is to be justified by 
the solution of Eq. ([2]) in the weak field limit, which in 
fact scales as 0, M 



The comparison to data obtained in the numerical simu- 
lation is presented in Fig.[T] This simple log-ansatz works 
well as long as g 2 /J,r < 0.5 and is stable under the a — > 
limit as shown by a thin curve in the figure. 

So far, we reached an ansatz (j 1 5[) with the infrared 
cut-off provided by the nucleus size L in a heuristic way. 
In reality, however, the long-ranged correlation should be 
cut off by the confining scale ~ Aq CD rather than L. In 
the continuum limit, hence, the initial energy density in 
the central collision at r < (g 2 /i) _1 should be estimated 
by 
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^( 5 V) 4 {ln(A Q ^ D /r)} . (16) 



In writing Eq. (|16p we put the constants c and tt ap- 
pearing in Eq. (|T5|l away into ambiguity of Aqcd- We 
remark that the AQCD-dependence would be milder than 
the above in the regime after the "formation time" as 
investigated in Ref. 

It is interesting to apply our formula (|16p to the Au-Au 
collision at RHIC in the physical unit. We make use of 
the parameter choice as commonly used in the numerical 
simulation, i.e., g 2 /An = 1/n and g 2 ^ = 2 GeV 0, H, @, 
[Toj . As for the confining scale, we vary Aq CD from 1 fm 
to 12 fm ~ L. The results are summarized as follows; 
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e(r = 0.1 fm) [ GeV fm" 3 ] 
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Our log-ansatz overestimates the energy density and the 
third row shows the corrected values with a factor 0.67 
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The energy density is dominated only by the Abelian 
part — {dA) 2 , hence it follows the tq/t scaling. We would 
comment on a curious observation that, if we extrapolate 
our estimate e(t = 0.1 fm) up to r = 1 fm assuming the 
tq/t scaling, the initial energy density obtained accord- 
ingly is very close to the standard estimate by means of 
the Bjorken formula; e(r = l fm) — 5.1 GeV- fm -3 . 

We shall next consider the problem of instability in 
the rest of this paper. We treat fluctuations SAi and 
8E l in the linear order around the boost-invariant CGC 
background which we discussed above. As formulated 
in Ref. [Hj], 5E n should be constrained by the Gauss 
law and we drop 5A V because it is accompanied by r 2 
from the metric. In what follows we regard ^-dependent 
fluctuations, 5Ai and 8E l , as the "soft" fields and rj- 
independent CGC fields as the "hard" background which 
brings about instability. The linearized equations of mo- 
tion are 



rd T 5A % = 8E l , 
d T S& = -t- 1 v 2 5A 1 



(18) 



where we introduced the Fourier transform 5Ai from r\ to 
the wave number v (i.e. d 2 5Ai(rj) — > — v 2 8Ai(v)) and we 
denoted the inverse of the transverse background gluon 
propagator as G" 1 whose definition is 
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We note that, in correspondence to the HEL approach, 
the color current encoding the anisotropic distribution 



4 




8 M T 

FIG. 2: Instability tendency with the initial CGC background 
fixed at r = in case of L/a = 700 for different wave numbers 
v = 1, 5, and 10 from the left to the right. 

of hard background is identified as jf = [G^ — 

5 ah (Sijdkdk — didj)~\ 5A°j. Although it is not clear whether 
this current could have anything to do with that in the 
HEL approach after the ensemble average, we can shortly 
confirm that instability may occur even in the purely clas- 
sical regime. 

It is easy to solve Eq. (fTgf to obtain SAi as a function 
of a given constant CGC background, at initial time, 

because we can diagonalize Gpl in a proper basis of 8A\. 
If we write its eigenvalues as A, the solution is 

6Ai(X) = cm ReIiv(V\T) + c 2l Im/ W (%/Ar) , (20) 

for A > (which exponentially grows) and 

SAi(X) = c Xi ReJi J/ ( v ^A|r) + c 2l ImJ^C^Ajr) , (21) 

for A < (which oscillatorily diminishes), where J n (x) 
and I n (x) are the first-kind and modified Bessel func- 
tions. These special functions are singular as (\/[AjT) ±w 
which furiously rotates in complex space as r — * 0. At 
later time when the asymptotic behavior (jTTJ) realizes 
due to expansion, A is no longer a constant but a func- 
tion of time like ~ £/r with some dimensionful constant 
£, which results in the solutions l2±iv(^V^r) for £ > 
and J2±w(2\/\£\t) for £ < 0. For other general cases, 
the r-dependence in the eigenvalue is intricate, and one 
needs to solve Eq. (jTSJ) numerically. 

We do not treat such a general case but what we will 
pursue here is to clarify whether the initial CGC fields 
could induce exponential growth for soft degrees of free- 
dom. We consider this because the initial CGC config- 
urations at r = are the most (spatially) anisotropic 
(namely, large longitudinal and zero transverse fields) 
and thus anticipated to cause the most unstable modes. 
Under such extreme circumstances we expect that the 



physics of instability becomes clear. Thus, as a first trial, 
it should be an appropriate starting point. We can say 
that what we will do is to extract the tendency toward in- 
stability under the presence of the CGC background. For 
simplicity we will focus on the case that fluctuations are 
uniform in transverse space, i.e. diSA = 0, which should 
be the most unstable, as adopted in Ref. [la |. 

We should calculate the Gaussian average ([5|) of the 
transverse chromo-magnetic field ((Bi) 2 ) ~ ((vSAi) 2 ) 
and the chromo-electric field ((Ei) 2 ) ~ ((rdrSAi) 2 ) 
which were as small as ~ r 4 previously but non-zero 
this time with fluctuations depending on n, contribut- 
ing to the longitudinal pressure The straightforward 
calculation is, however, technically hard. We will approx- 
imately do it by picking the mean value up, 

<S)« ) = 5 ^ abl = -M" 6 ! 47 (2 2 M) 2 , (22) 

where a is defined in Eq. |7]). and taking the ensemble 
average over its dispersion, 

(23) 

where x is defined in Eq. (| 1 2[) . 

Because A is negative, the soft fluctuations in the vicin- 
ity of the averaged CGC background are stable belong- 
ing to the type of solution (f2~Tj) . The eigenvalue of G^ 
distributes according to random CGC configurations and 
spreads from A with the dispersion <5A, meaning that some 
CGC configurations may have negative A. That is, if we 
evaluate, 

/oo 
dXO[8A(X)]e-^-^' 2SX \ (24) 
-oo 

using Eq. (|20p . the contributions near A ~ A dominate 
only when time is small until the negative A constituents 
grow up as time elapses. The transverse field strengths 
obtained in this way are plotted in Fig. [2] 

To draw Fig. [2J we chose the initial time g 2 /J.To = 0.001 
at which we set c\ and c 2 of Eq. (|20[) or (f2Tj) by the ini- 
tial condition, 5Ai = cjy/v and SEi = Tod T 5Ai = C\/v, 
inspired by quantum fluctuations discussed in Ref. [lq |. 
It is interesting to see that this specific initial condition 
(SEi - makes ((Bi) 2 ) ~ ((vSAi) 2 ) and ((Ei) 2 ) 

comparable to each other, leading to their almost alter- 
nate oscillations. That is why the sum of transverse field 
strengths depicted in Fig. [2] never come close to zero, 
which makes a contrast to the results in Ref. [151 ]. 

We can conclude that there is certainly the tendency 
toward instability associated with initial CGC back- 
ground. The onset of instability in the present case is 
located much earlier than preceding works. It is be- 
cause we only investigated the strongest instability en- 
compassed in the initial CGC fields. Because the CGC 
background itself evolves with time, in fact, the genuine 
growth of instability should be slower and weaker than 
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shown in Fig. [2l Nevertheless, we can learn the qualita- 
tive character of the "Glasma" instability. The intuitive 
picture is as follows. The soft fluctuations of gluon fields 
are non-Abelian charged and feel a force under influence 
from the CGC background. The ensemble of random 
CGC distribution contains not only color fields which 
suppress the color current provided by charged soft fluc- 
tuations but also color fields which amplify the current. 
Although the current is suppressed on average, the large 
r behavior is predominantly determined by mixture of 
CGC fields which enhance the input current. Therefore, 
we think that it is rare fluctuation in the CGC ensemble 
from which the Glasma instability can occur. 

It is necessary to deal with A(r) as not a constant but 
a function of r in order to quantify instability further. 
In our treatment mentioned above we dropped the effect 
of longitudinal expansion for the hard part, while the ex- 
ponential growth should take a form of l2±iv{1\f^r) ~ 
r -1 / 2 cxp[2 v / ^r] asymptotically when A ~ £/t, as we re- 
marked before. The analytical estimate of £ deserves 
future clarification. Also, we have to evaluate A and SX 
in a resummed form like Eq. (|15[) beyond the naive ex- 
pressions (j2"2")l and (f2U)) . As a matter of fact, the growth 
rate seems to be determined by A and 5 A regardless of v 
in view of our results in Fig. [2] Quantitative details of an 
analytical description should be improved with guided by 
systematic instability studies in the numerical simulation 



in the future. 

In summary, we developed an analytic formula to es- 
timate the initial energy density. Our conclusion is 
e(r = 0.1 fm) = 40 - 50 GeV- fm~ 3 in the (central) Au- 
Au collision at RHIC. The uncertainty comes from the 
infrared cut-off (or confining) scale. Also, we analyzed 
the tendency toward instability in the presence of the 
initial CGC background fixed at r = 0. We found that 
there exist unstable modes as a result of the ensemble 
average of random CGC configurations, some of which 
strengthen the color current brought in by soft fluctua- 
tions. Although the Glasma instability might have a con- 
nection to non-Abelian plasma instabilities at a deeper 
level, we would emphasize that we could understand the 
Glasma instability not relying on the picture of plasma 
instabilities that are premised on anisotropic momentum 
distribution of hard particles. The bottom line is, thus, 
that the Glasma instability exists from r = even when 
the particle picture is irrelevant yet. 
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